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I . SUMMARY 

The  objective  of  our  research  program  is  to  exaunine 
the  parameters  that  affect  the  seismic  signals  from  under- 
ground explosions.  Our  attention  is  primarily  directed  to 
those  features  of  the  seismic  waveforms  that  reliably  indi- 
cate the  explosion  yield.  Our  research  program  includes 
empirical  studies  of  the  available  data,  experimental  studies 
using  small  charges  to  simulate  explosions  and  the  develop- 
ment and  application  of  theoretical  and  numerical  methods. 
Emphasis  is  on  the  latter,  in  particular  we  are  applying 
techniques  for  numerically  simulating  the  far-field  signals 
from  both  contained  and  cratering  underground  explosions. 

The  numerical  simulation  techniques  represent  a synthesis 
of  the  finite  difference  methods  for  computing  ground'motion 
in  the  near-source  large  displacement  regime  and  the  effi- 
cient wave  propagation  techniques  of  theoretical  seismology. 

During  the  first  three-month  period  of  our  present 
contract,  our  research  activities  have  been  concentrated  in 
five  areas.  These  areas  and  the  major  results  obtained  in 
each  are  summarized  in  the  next  five  paragraphs. 

We  have  been  involved  in  a joint  project  with  Applied 
Theory,  Inc.  (ATI)  and  Pacific  Sierra  Research  (PSR)  in 
which  the  main  objective  has  been  to  study  the  variation  of 
the  teleseismic  magnitudes,  M^  and  m^^,  with  source  material 
and  depth  for  a series  of  cratering  calculations.  Our  con- 
tribution has  been  to  process  the  output  of  the  ATI  calcula- 
tions and  obtain  estimates  of  M^  and  rtij^.  A detailed  report 
is  being  prepared  describing  our  techniques  and  results. 

In  Section  3.1  of  this  report  we  present  the  synthetic 
seismograuns  and  the  and  M^  values  we  have  computed.  These 
results  are  presented  without  comment  regarding  their  signi- 
ficance or  importance  for  treaty  monitoring.  Such  comments 
will  be  more  appropriate  after  detailed  consultation  with 
ATI  and  PSR  scientists  and  will  appear  in  our  topical  report. 
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in  important  question  now  under  study  at  many  institu- 
tions is  the  possibility  that  NTS  may  be  underlain  by  a re- 
gion of  anomalously  high  attenuation  which  causes  the  m^^  from 
NTS  explosions  to  be  biased  low.  We  have  also  been  studying 
this  question.  As  a first  step  we  have  compared  the  frequency 
content  of  a large  population  of  NTS  and  Eurasian  explosions 
at  several  recording  stations  where  we  have  previous  experience 
in  synthesizing  short  period  recordings  of  NTS  events.  Our 
results  are  shown  graphically  and  discussed  in  Section  3.2. 

We  conclude  that  the  dominant  periods  for  Kazakhstan  events 
are  considerably  shorter  than  for  NTS  events  in  the  saune  yield 
rcinge.  The  differences  are  as  much  as  0.5  seconds  at  two  of 
the  stations.  A significant  point  about  the  NTS  data  is  that 
the  period  of  the  PILEDRIVER  event,  an  unusual  NTS  explosion 
in  granite,  falls  within  the  NTS  population.  While  these 
data  do  not  argue  conclusively  for  differences  in  attenuation 
being  the  cause  of  the  observed  frequency  differences,  they 
point  in  that  direction,  particularly  in  view  of  the  large 
sampling  of  event  depths  and  source  media.  We  are  continuing 
to  study  this  problem  using  our  synthetic  seismogreun  generat- 
ing capabilities. 

r 

k-  ! An  important  part  of  our  research  capability  is  our 

numerical  (finite  difference)  techniques  for  simulating  the 
coupling  of  the  explosion  energy  into  elastic  waves  in  the 
earth.  Using  our  one-dimensional  (spherically  symmetric) 
computer  code,  we  have  been  studying  the  effect  of  near 
source  material  properties  on  the  seismic  signal.  Our  pre- 
viously developed  constitutive  models  were  found  to  be  un- 
satisfactory when  used  for  highly  porous  materials.  Our 
constitutive  models  were  reevaluated  and  the  results  are  dis- 
cussed in  Section  3.3.  A minor  modification  to  the  consti- 
tutive model  was  found  to  enable  the  program  to  give  ex- 
cellent results  for  highly  porous  materials.  The  new  model 
is  as  compatible  with  laboratory  data  as  was  the  one  pre- 
viously used. 


2 


R-3108 


-■T* 


During  the  quarter,  a theoretical  study  was  conducted 
in  which  the  objective  was  to  compare  explosions  in  salt  in 
Eurasia  to  similar,  though  hypothetical,  events  at  NTS  and 
to  granite  events  at  NTS.  The  results  of  this  study  were 
presented  in  detail  in  a topical  report  by  Bache,  Cherry 
and  Mason  [1976] . The  main  results  from  that  report  are 
given  here  in  Section  3.4.  A series  of  equivalent  point 
source  representations  for  explosions  at  various  depths  in 
salt  were  computed.  We  found  that  the  long  period  level, 
proportional  to  M^,  was  only  mildly  sensitive  to  depth.  The 
short  period  level,  however,  was  much  more  sensitive,  varying 
by  a factor  of  two  for  depths  between  150  and  2500  meters. 

This  is  due  to  the  opening  of  tensile  cracks  in  the  material 
at  the  low  confining  pressures  characteristic  of  shallow 
depths.  We  found  that  when  all  other  factors  are  the  same, 
explosions  in  salt  couple  about  a factor  of  two  smaller  than 
events  in  PILEDRIVER-type  granite.  We  also  looked  at  the 
effect  of  attenuation  on  the  seismograms.  If  we  assume  the 
attenuation  for  NTS  events  is  about  40  percent  higher  than 
for  Eurasian  events,  we  find  that  identical  salt  explosions 
at  the  two  test  sites  would  give  body  wave  amplitudes  that 
differ  by  a factor  of  = 3-4. 

We  have  been  engaged  in  an  experimental  program  in 
which  small-scale  modeling  experiments  of  underground  explo- 
sions are  to  be  compared  to  computer  simulations  of  the 
experiments.  After  some  difficulties  with  the  explosives 
and  displacement  gauges  which  took  several  months  to  resolve, 
several  successful  modeling  experiments  have  been  carried 
out.  The  test  medium  is  a concrete  cylinder  and  the  charges 
are  0.25  grams  of  PETN  explosive  which  can  be  buried  at 
selected  depths.  In  Section  3.5  we  give  some  results  from 
our  first  experiments.  These  results  indicate  the  reliability 
and  repeatability  of  the  experiments  and  are  necessary  to 
lend  credibility  to  the  results.  The  ground  motion  records 
are  of  excellent  quality  and  we  are  now  proceeding  to 
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numerically  model  these  experiments.  These  should  provide 
an  excellent  normalization  of  our  theoretical  techniques. 
Further  experiments  are  being  performed  in  which  the  charge 
will  be  shallow  so  cratering  will  result. 
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II.  INTRODUCTION 

The  primary  objective  of  the  research  prograun  is  to 
systematically  examine  the  parameters  that  control  magnitude- 
yield  relationships  for  underground  nuclear  explosions.  During 
the  first  three  month  period  of  this  contract,  activity  has 
been  concentrated  in  the  following  areas: 

1.  The  calculation  of  teieseismic  magnitudes  (M^ 
and  mj^)  for  nuclear  explosion  cratering  calcula- 
tions carried  out  by  Applied  Theory,  Incorporated. 

2.  The  study  of  body  wave  attenuation  properties  for 
selected  non-U. S.  travel  paths. 

3.  Development  of  improved  techniques  for  one- 
dimensional numerical  simulation  of  explosions. 

4.  A study  of  the  dependence  of  body  wave  magnitude 
on  yield  for  underground  explosions  in  salt. 

5.  A laboratory  simulation  of  underground  explosions 
(project  carried  over  from  our  previous  contract. 
Project  VT/6712). 

The  following  section,  the  Technical  Discussion,  is  de- 
voted to  technical  summaries  of  our  work  in  each  of  the  above 
listed  areas. 
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III.  TECHNICAL  DISCUSSION 

3.1  CALCULATION  OF  FOR  ATI  CRATERING  CALCULATIONS 

3.1.1  Introduction 

We  are  involved  in  a joint  project  with  Applied  Theory, 
Incorporated  (ATI)  and  Pacific  Sierra  Research  (PSR)  to  study 
the  variation  of  the  teleseismic  magnitudes,  M and  m,  , with 

S iD 

source  material  and  depth  for  a series  of  cratering  calcula- 
tions. The  S*  function  has  been  to  process  the  output  of  the 
finite  difference  calculations  performed  by  ATI  and  to  deter- 
mine Mg  and  m^^  for  each  calculation.  This  is  a complex  and 
rather  difficult  task  and  a detailed  report  is  being  prepared 
describing  our  techniques  and  the  results.  In  this  report  we 
will  present  the  most  important  results.  We  would  li’.e  to 
point  out  that  the  results  are  preliminary  and  some  chan>  es 
may  result  after  consultation  with  ATI  and  PSR  scientists. 

3.1.2  Computational  Procedure 

Briefly  outlined,  the  computational  procedure  followed 
to  obtain  and  for  the  cratering  calculations  is  as  fol- 
lows : 

1.  The  ATI  performed  cratering  calculations  were  carried 
into  the  regime  where  the  material  response  is  approxi- 
mately linearly  elastic.  A tape  was  then  prepared 
containing  the  time  histories  of  che  divergence  and 
curl  (7*u(t),  7xu(t))  of  the  displacement  field  on  a 
radius,  denoted  the  elastic  radius,  centered  at  ground 
zero.  The  geometry  and  coordinate  system  are  shown 

in  Figure  3.1.1. 

2.  The  divergence  and  curl  are  expanded  in  a series  of 
spherical  harmonics  to  obtain  an  equivalent  elastic 
source.  The  procedure  is  formally  that  described  by 
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Bache  and  Harkrider  [1976]  for  sources  in  a whole 
space.  However,  the  presence  of  a free  surface  re- 
quires a number  of  assumptions  that,  to  some  extent, 
control  the  solution.  The  extent  of  this  effect  is 
discussed  in  some  detail  in  our  forthcoming  report. 
However,  the  most  important  point  is  that  all  the 
cratering  calculations  are  treated  the  same  way  and 
the  relative  values  of  tiij^  and  should  be  preserved. 

Using  the  equivalent  elastic  source,  synthetic  seis- 
mograms are  computed  for  body  waves  and  surface 
waves.  For  body  waves  the  pertinent  references  are 
Bache  and  Harkrider  [1976]  and  Bache,  et  al.  [1976]. 
For  surface  waves  we  use  the  method  of  Harkrider 
which  has  been  described  in  numerous  publications; 
e.g.,  Harkrider  [1964],  with  certain  modifications 
indicated  in  Section  3.1.6. 

The  values  of  itij^  and  are  strongly  dependent  on 
the  crust  and  upper  mantle  models  used  in  the  syn- 
thetic seismogram  calculations.  However,  the  same 
models  are  used  in  all  calculations  and  the  relative 
values  should  be  insensitive  to  these  models. 

The  free  surface  stresses  due  to  the  ejecta  fallback 
were  also  provided  to  us  by  ATI.  These  data  were 
analyzed  to  determine  the  effect  on  the  far-field 
body  and  surface  waves. 


j.i.  j Description  of  Cratering  Calculations 

We  will  be  describing  fourteen  source  calculations  in 
three  emplacenent  materials.  Two  of  these  are  one-dimensional 
calculations  for  contained  explosions  in  a homogeneous  whole 
space  and  provide  benchmark  cases  for  measuring  the  effect  of 
the  cratering.  The  other  twelve  are  cratering  calculations. 
The  important  parameters  characterizing  the  calculations  are 
summarized  in  Table  3.1.1. 


I 
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TABLE  3.1.1 

PARAMETERS  DESCRIBING  THE  ATI  CALCULATIONS 


Identifier 

Material 

D^Jth  (km) 

aCknv/sec) 

6(kiTv/sec) 

p(kiiv'sec) 

1 

Granite 

0.159 

4.402 

2.54 

2.661 

2 

n 

0.207 

4.402 

2.54 

2.661 

3 

n 

0.253 

4.406 

2.542 

2.661 

4 

Dry  Semdstone 

0.159 

2.822 

1.740 

2.30 

5 

•f 

0.207 

2.825 

1.743 

2.30 

6 

tt 

0.253 

2.828 

1.744 

2.30 

7 

11 

0.207 

2.836 

1.755 

2.30 

8 

Wet  Sandstone 

0.053 

2.620 

1.509 

2.40 

9 

n 

0.159 

2.624 

1.592 

2.40 

10 

n 

0.207 

2.614 

1.513 

2.40 

11 

n 

0.253 

2.618 

1.517 

2.40 

12 

It 

0.531 

2.619 

1.519 

2.40 

13 

Granite 

Spherically 

Synmetric 

4.239 

2.448 

2.661 

14 

Wet  Sandstoie 

SphericcC-ly 

Syrtinetric 

2.619 

1.530 

2.40 

9 
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We  see  from  the  table  that  the  calculations  have  the 
potential  to  help  us  understand: 

1.  The  effect  of  burial  depth  on  the  teleseismic 
signature  of  cratering  explosions. 

2.  The  effect  of  emplacement  material  on  the  signal 
from  cratering  explosions. 

3.  The  differences  between  cratering  and  contained 
shots  in  the  same  material. 

Considering  the  approximations  made  in  our  calculations  (ex- 
clusive of  any  difficulties  with  the  finite  difference  cal- 
culations) , the  above  are  listed  in  order  of  the  amount  of 
confidence  we  place  in  the  results. 


s' 


3.1.4  Far-field  Displacement  Spectra 

For  each  of  the  fourteen  sources  listed  in  Table  3.1.1 
we  compute  far-field  displacement  spectra.  For  the  one- 
dimensional, spherically  symmetric  calculations  the  equivalent 
elastic  source  is  the  reduced  displacement  potential,  'i' (t-R/a)  , 
which  is  related  to  the  displacement  by 


u (R,a)) 


(3.1.1) 


where  all  quantities  are  F irier  transformed.  Then  the  far- 
field  displacement  spectr  .s  defined  by 


Ujj  (R,u)) 


(3.1.2) 


An  entirely  analogous  procedure  is  followed  for  the  com- 
plex equivalent  elastic  source  representations  of  the  cratering 
calculations.  That  is,  we  retain  only  terms  of  order  R~  in 
the  expansion. 
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The  displacement  spectra  are  presented  in  Figures 
3.1.2  - 3.1.7.  The  plots  are  log-log  in  cunplitude  versus 
frequency.  Note  that  the  scale  on  the  amplitude  axis  is  in 
powers  of  10  while  the  actual  frequencies  are  printed  on 
the  abscissa. 

Two  frequencies  are  singled  out  on  the  spectra  for 
the  cratering  calculations  and  are  denoted  A and  B.  The 
frequency  denoted  B is  associated  with  the  (approximate) 
total  time  of  the  ground  motion  data  provided  by  ATI.  That 
is,  if  the  V*u(t)  and  Vxu(t)  at  a typical  station  on  the 
elastic  radius  had  their  first  non-zero  value  at  t^  = 0.4 
seconds  and  the  final  time  point  was  at  t^  = 2.5  seconds,  we 
say  that  the  frequency  B is  l/Ct^-t^^)  or,  for  this  example, 
0.48  Hz.  Then  B is  the  lowest  frequency  that  can  be  assoc- 
iated with  the  actual  computed  data. 

The  meaning  of  the  frequency  A is  more  difficult  to 
explain  and  the  interested  reader  is  referred  to  our  forth- 
coming detailed  report.  Beyond  the  time  t^  we  assume  that 
the  lowest  order  terms  (the  dipole)  in  the  expansion  of  the 
divergence  and  curl  are  zero.  The  higher  order  terms  are 
assumed  to  remain  static  at  the  value  reached  at  the  last 
time  point.  The  amplitudes  at  all  frequencies  below  B are 
dependent  on  these  assumptions.  In  order  to  compute  the 
Fourier  transform  we  extend  the  time  histories  out  to  some 
t^  using  the  assumptions  mentioned  above.  Then  A a l/(tg-t^). 

How  do  we  compute  values  for  frequencies  below  A?  We 
know  that  the  response  at  low  frequencies  is  dominated  by  the 
dipole  term  in  the  expansion.  Further,  we  can  prove  that  the 
dipole  term  must  behave  like  u*  at  low  frequency  for  a bounded 
solution.  This  corresponds  to  a flat  far-field  displacement 
spectrum  at  low  frequency.  Our  actual  procedure  is  to  extra- 
polate the  amplitude  of  the  dipole  term  from  its  value  at  A 
by  assuming  proportionality  to  u*.  This  is  nearly  the  same 
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Figure  3.1.4.  Far-field  displacement  spectra  for  four  cratering  explosions  In  dry 
sandstone.  The  spectra  are  for  a takeoff  angle  t = 20*. 
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as  assuming  the  far-field  displacement  spectrum  to  be  flat 
to  long  periods  from  its  value  at  A as  can  be  seen  in  the 
figures . 


Finally,  in  Tables  3.1.2  - 3.1.3  we  tabulate  the 
spectral  values  that  have  the  greatest  significance  for  com- 
puting teleseismic  body  and  surface  waves.  The  amplitude 
values  in  Tables  3.1.2  - 3.1.3  are  all  on  the  same  scale 

. 4 

(all  have  been  multiplied  by  10  R)  and  are  best  viewed  as 
relative  amplitudes.  The  reliability  of  these  values  for 
scaling  itij^  and  is  discussed  in  our  topical  report. 


3.1.5  Body  Wave  Amplitude,  m^^ 

In  this  section  we  present  our  synthetic  seismograms 
for  the  fourteen  ATI  sources  and  give  the  values  for  each. 
The  synthetic  seismogram  calculations  include  the  following 
elements : 

1.  The  equivalent  elastic  sources  which  give  the  (whole 
space)  far-field  displacement  spectra  described  in  the 
previous  section  are  embedded  in  a layered  model  of 
the  crust  in  the  source  region.  The  basic  model  used 
for  the  calculations  is  tabulated  in  Table  3.1.4. 

The  top  layer  was  changed  to  have  the  properties  ap- 
propriate to  each  source  as  listed  in  Table  3.1.1. 

In  carrying  out  the  calculations  only  the  downgoing 
waves  emitted  by  the  source  are  computed;  that  is, 
no  free  surface  is  included  in  the  source  crustal 
model. 

2.  The  far-field  body  waves  emanating  from  the  base  of 
the  source  crust  and  characterized  by  ray  parameter 
p = 0.079  sec/km  are  calculated. 

3.  The  upper  mantle  is  accounted  for  by  a step  function 
response  computed  using  generalized  ray  theory.  In 
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TABLE  3.1.2 

BODY  WAVE  AMPLITUDE  SPECTRA  - FAR-FIELD  P WAVE 
1 HZ  SPECTRAL  AMPLITUDES  AT  T = 20“ 


Identifier  Material 


Depth  (km) 


Spectral  Amplitude 


13 

Granite 

Spherically 

Symmetric 

00 

• 

1 

n 

0.159 

1.8 

2 

n 

0.207 

2.1 

3 

n 

0.253 

1.5 

4 

Dry  Sandstone 

0.159 

0.60 

5 

n 

0.207 

0.65 

6 

ti 

0.253 

0.72 

7 

Weak  Dry  Sand- 
stone 

0.207 

1.4 

14 

Wet  Sandstone 

Spherically 

Symmetric 

1.7 

8 

n 

0.053 

1.3 

9 

n 

0.159 

2.1 

10 

It 

0.207 

2.8 

11 

It 

0.253 

2.7 

12 

II 

0.531 

2.5 

19 


TABLE  3.1.3 

SURFACE  WAVE  AMPLITUDE  SPECTRA  (0.05  Hz) 


Identifier 

P (t=20‘’) 

S(t=20“) 

P(t=70‘’) 

S (t=70 

13 

1.3 

- 

1.3 

- 

1 

2.9 

2.0 

1.1 

5.4 

2 

3.2 

5.7 

1.2 

15.6 

3 

3.3 

0.9 

1.2 

2.7 

4 

6.5 

1.1 

2.4 

2.8 

5 

7.1 

5.0 

2.6 

1.4 

6 

5.6 

5.7 

2.1 

1.6 

7 

10.6 

0.95 

3.9 

2.4 

14 

2.0 

- 

2.0 

- 

8 

2.2 

1.7 

0.78 

4.6 

9 

10.1 

7.4 

3.7 

20.2 

10 

21.3 

13.6 

7.7 

37.8 

11 

14.4 

22.0 

5.2 

60.1 

12 

13.7 

11.8 

5.0 

32.5 

20 
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TABLE  3.1.4 

SOURCE  REGION  CRUSTAL  MODEL 


Depth  (kn) 

Thickness  (km) 

a (km/sec) 

6 (km/sec) 

P (g/cm 

1.0 

1.0 

Granite, 

Wet  or  Dry  Sandstone 

1.7 

0.7 

4.7 

2.7 

2.6 

to 

• 

1.0 

5.4 

2.8 

2.7 

4.0 

1.3 

5.8 

3.45 

2.8 

20.0 

16.0 

6.0 

3.50 

2.8 

TABLE  3.1.5 

RECEIVER  REGION  CRUSTAL  STRUCTURE 


Depth  (kn) 

Thickness  (km) 

a ()cm/sec) 

6 (km/sec) 

p (km/sec) 

2.58 

2.58 

3.67 

2.31 

2.40 

4.84 

2.26 

5.42 

3.27 

2.60 

11.61 

6.77 

5.80 

3.45 

2.60 

20.0 

8.39 

6.00 

3.50 

2.80 

21 


this  case  we  took  the  distance  to  be  A = 36“  which 
is  beyond  the  upper  mantle  triplications  and  the 
upper  mantle  response  is  essentially  a constant 
geometric  spreading  factor. 

4.  The  response  of  the  receiver  crustal  model  (Table 
3.1.5)  is  included.  In  this  case  the  receiver 
crust  has  little  effect  other  than  scaling  the 
seismogram. 

5.  The  response  is  convolved  with  an  operator  repre- 
senting the  attenuation  and  dispersive  properties 
of  the  earth.  For  these  calculations  we  took 

t*  = 0.7  (t*  is  the  ratio  of  travel  time  to  the 
effective  path  attenuation  factor,  Q) . 

6.  The  ground  motion  is  convolved  with  the  response 
of  a standard  short  period  seismograph  system. 


The  synthetic  seismograms  are  shown  in  Figures  3.1.8  - 
3.1.10.  The  cycle  from  which  is  measured  is  indicated  on 
each  record  by  a bar.  The  itij^  values  and  the  period  of  the 
cycle  are  tabulated  in  Table  3.1.6.  The  nij^  is  computed  from  |' 

mj^  = log  I + 3.32  , (3.1.3) 

where  T is  the  period  tabulated,  A is  the  peak-to-peak  ampli- 
tude of  the  indicated  cycle  corrected  for  the  instrument 
response  at  the  period  T,  and  the  constant  3.32  is  the  ap- 
propriate distance  correction  factor.  Recall  that  for  all 
the  seismogram  calculations  there  is  no  free  surface  near 
the  source.  i< 

It  is  interesting  to  examine  how  the  itij^  scales  with 


the  far-field  displacement  spectra  discussed  in  Section 
To  first  order  we  expect  to  have  (Bache,  et  al. , 1976] , 

3.1.4. 

mb  = log  (“s  “)  ' 

(3.1.4) 

j 
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Figure  3.1.8. 


TIME  (SEC) 

Synthetic  short  period  seismograms  for  one 
spherically  symmetric  and  three  cratering 
calculations  in  granite.  The  numbers  to  the 
left  are  ground  motion  in  microns  at  1 Hz. 


Identifier 


,202*00 


.273*00 
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12<J'*00 

,244*00 


,251*00 

.344*00 


.333*00 

,223*00 


,223*00 
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Figure  3.1.10. 


TIME  (SEC) 

Synthetic  short  period  seismograms  for  one 
spherically  symmetric  and  five  cratering 
calculations  in  wet  sandstone.  The  numbers 
to  the  left  are  ground  motion  in  microns  at 
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TABLE 

3.1.6 

SUMI4ARY 

OF  Itij^  VALUES 

FOR  ATI  CRATERING 

CALCULATIONS 

Identifier 

Depth 

Period 

Normalized 

(A/D/OgU) 

Granite 

13 

Spherically 

Syninetric 

6.29 

0.81 

0.93 

1 

0.159 

6.28 

0.78 

0.85 

2 

0.207 

6.34 

0.82 

0.73 

3 

0.253 

6.28 

0.80 

1.01 

Dry  Sandstone 

4 

0.159 

5.48 

0.96 

0.96 

5 

0.207 

5.57 

0.88 

1.09 

6 

0.253 

5.55 

0.86 

0.94 

7 (weak) 

0.207 

5.80 

0.93 

0.86 

Wet  Seindstone 

14  (1-D) 

Spherically 

Synmetric 

5.88 

0.79 

1.00 

8 

0.053 

5.66 

0.94 

0.79 

9 

0.159 

5.98 

01^2 

1.02 

10 

0.207 

6.11 

0.92 

1.05 

11 

0.253 

6.10 

1.03 

1.05 

12 

0.531 

6.15 

0.99 

1.27 

26 
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where  u is  the  displacement  spectrtim  at  the  controlling  fre- 
quency and  a is  the  P wave  velocity  at  the  source.  Using 
s 

the  values  of  a from  Table  3.1.1  and  u (1  Hz)  from  Table 
s 

3.1.2,  we  normalize  the  amplitudes  and  scale  to  the  value 
for  calculation  14.  The  normalized  scaled  eunplitudes  are 
tabulated  in  the  last  column  of  the  table.  We  see  that  all 
values  are  within  + 27  percent  of  the  normalized  amplitude 
for  calculation  14.  Further,  eleven  of  the  fourteen  are 
scaled  to  within  15  percent  by  (3.1.4).  The  discrepancies 
are  mainly  attributed  to  the  fact  that  we  are  normalizing 
to  the  spectral  amplitude  at  1 Hz  while  the  frequencies  con- 
trolling the  seismogram  amplitudes  range  from  0.97  - 1.32  Hz. 
From  Figures  3.1.2  - 3.1.7  we  see  that  the  P wave  spectra  are 
complex  and  rapidly  changing  in  this  region. 

3.1.6  Surface  Wave  Amplitude, 

In  this  section  we  present  our  synthetic  long  period 
seismograms  for  the  fourteen  ATI  sources  and  give  the 
values  for  each.  The  synthetic  seismogram  calculations 
include  the  following  elements; 

1.  The  method  for  the  surface  wave  calculations  is 
described  by  HarJcrider  [1964]  and  Har)crider  and 
Archambeau  [1977] . The  same  equivalent  elastic 
source  formulation  used  for  the  body  waves  is  used 
for  the  surface  waves.  Once  again,  only  the  down- 
going waves  from  the  source  are  included  in  the 
calculations. 

2.  Two  crustal  models  are  used  for  the  path,  one  for 
the  very  near  source  region  and  one  for  the  re- 
mainder of  the  path  to  the  receiver.  The  average 
path  model  is  one  proposed  for  North  America  by 
McEvilly  [1964] . The  only  difference  between  the 
models  is  that  the  top  three  kilometers  of  the 
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source  region  crustal  model  is  replaced  by  the  ATI 
granite,  wet  or  dry  sandstone,  depending  on  the  source 
material.  For  the  long  periods  controlling  teleseismic 
Mg  the  reflection  coefficient  for  Rayleigh  waves 
passing  across  the  boundary  between  the  source  and 
average  path  crustal  model  is  close  to  unity. 

3.  The  ground  motion  is  convolved  with  the  response  of 
an  LRSM  long  period  seismometer.  A Q operator  which 
has  only  a minor  effect  is  also  included.  The  Q 
model  is  that  of  Tryggvason  [1965]  . 

4.  The  seismograms  were  synthesized  at  a range  of 
3000  km.  The  was  computed  using  the  formulas  of 
Marshall  and  Basham  [1972] . For  this  range  the 
formula  reduces  to 

M = log  A + 1.38  + P(T) 
s 

where  A is  the  maximum  amplitude  (zero-to-peak)  of 
the  signal  with  period  near  20  seconds  and  P(T)  is 
a period  dependent  path  correction  tabulated  by 
Marshall  and  Basham  [1972].  The  correction  is  quite 
small  for  periods  near  20  seconds. 

The  vertical  component  Rayleigh  wave  synthetic  seismo- 
grams are  shown  in  Figures  3.1.11  - 3.1.13.  Two  seismograms 
are  shown  for  the  spherically  symmetric  contained  explosions 
in  granite  and  wet  sandstone  (cases  13  and  14).  In  one  of 
these  (13b,  14b)  only  the  downward  waves  are  included  in  the 
calculation.  The  seismograms  are  then  the  analog  of  the 
case  13  and  14  body  wave  seismograms  of  Figures  3.1.8  and 
3.1.10.  Since  the  M is  relatively  insensitive  to  depth 

S 

for  contained  explosions  (unlike  which  is  strongly  af- 
fected by  the  pP  phase) , we  also  computed  seismograms  for 
fully  contained  explosions  at  a depth  of  200  meters.  The 
seismograms  are  givrn  as  cases  13a  and  14a. 
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47?*01 


Identifier 


13a 

.■-.yVvWA 

(total  source)  f 1/ 


440-*01 


210*01 


1?4‘01 

127*01 


(down'^ar^  ^ V 

waves  only) 


171*01 


.S57*01 

.rf2?*01 


.525*01 


0.  50.  100.  150.  200.  250.  300.  350.  400.  450.  500.  550. 

TIME  (SEC) 


Figure  3.1.11. 


Vertical  component  Rayleigh  wave  seismograms 
for  a spherically  symmetric  and  three  cratering 
explosions  in  granite.  The  numbers  at  the  left 
are  displacement  in  microns  at  25  seconds. 


■.57tf*01 


0.  50.  100.  150.  200.  250.  300.  350.  400.  450.  500.  550 

TIME  (SEC) 

Figure  3.1.12.  Vertical  component  Rayleigh  wave  seismograir 


Vertical  component  Rayleigh  wave  seismograms 
for  four  cratering  explosions  in  dry  sand- 
stone. The  numbers  at  the  left  are  dis- 
placement in  microns  at  25  seconds. 
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Identifier 


214*01 


202*01 

17X*01 


(downward 
waves  only) 


214*01 

151*01 


A ,A  r A 


147*01 
51  ^>*01 


5<J1*01 

10P*02 


122*02 


-.132*02 

.133*02 
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9 

— f'' 

4.71 

N-  \>  V/  V \J' 

i 

r 

5.00 


-.104*02  1 Li • 1 ♦ 

0.  50.  100.  150.  200.  250.  300.  350.  400.  450.  500.  550. 


Figure  3.1.13, 


TIME  (SEC) 

Vertical  component  Rayleigh  wave  seismograms 
for  a spherically  symmetric  and  five  cratering 
explosions  in  wet  sandstone.  The  numbers  at 
the  left  are  displacement  in  microns  at  25 
seconds. 
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TABLE  3.1.7 

SUMMARY  OF  M„  VALUES  FOR  ATI  CRATERING  CALCULATIONS 
s 


Identifier 

Depth  (km) 

M 

s 

M^  Period 
s 

Log  Spectral 
Anpljitude, 
log  A (20  sec) 

Mg-log 

Granite 

13A 

1-D  Total  Source 

4.57 

20.1 

1.92 

2.65 

13B 

1-D  IDown  Waves 
Only 

4.34 

21.5 

1.62 

2.72 

1 

0.159 

4.01 

20.0 

1.30 

2.71 

2 

0.207 

4.95 

23.6 

2.25 

2.70 

3 

0.253 

4.72 

21.8 

2.03 

2.69 

Dry  Sandstone 

4 

0.159 

4.66 

22.0 

1.92 

2.74 

5 

0.207 

4.68 

22.7 

1.92 

2.76 

6 

0.253 

4.60 

22.5 

1.83 

2.77 

7 (weak) 

0.207 

4.83 

22.0 

2.07 

2.76 

Wet  Scmdstone 

14A 

1-D  Total  Source 

4.17 

22.8 

1.52 

2.65 

14B 

1-D  Dcwn  Waves 
Only 

4.26 

20.0 

1.53 

2.73 

8 

0.053 

4.19 

21.2 

1.43 

2.76 

9 

0.159 

4.71 

21.2 

2.02 

2.69 

10 

0.207 

4.99 

19.3 

2.35 

2.64 

11 

0.253 

5.07 

20.1 

2.44 

2.63 

12 

0.531 

5.00 

19.8 

2.34 

2.66 

32 
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The  important  data  from  the  seismograms  of  Figures 
3.1.11  - 3.1.13  are  summarized  in  Table  3.1.7.  The  phase 
at  which  the  eimplitude  for  was  measured  is  indicated  on 
the  seismograms  of  Figures  3.1.11  - 3.1.13  by  a bar.  The 
period  of  this  phase  is  given  in  the  table.  Also  listed  in 
the  table  is  the  spectral  amplitude  of  the  true  ground  motion 
at  a period  of  20  seconds.  The  difference  between  and 
log  A indicates  the  consistency  to  which  the  measurement 
represents  a true  measurement  of  the  energy  at  frequencies 
in  this  range.  These  differences  are  tabulated  in  the  last 
column  of  Table  3.1.7.  We  see  that  the  values  are  quite 
consistent  with  the  spectral  measurements  with  the  spread 
between  the  maximum  and  minimum  values  being  0.14  units. 

3.1.7  Contribution  of  the  Ejecta  Fallback  to  Teleseismic 

Ground  Motion 

In  carrying  out  the  cratering  calculations,  ATI  also 
tracked  the  ejecta  and  computed  the  surface  loading  due  to 
its  fall  back  to  the  free  surface.  The  resulting  surface 
loading  consists  of  an  axisymmetric  distribution  of  normal 
and  shear  stresses  on  the  surface.  If  we  Fourier  transform 
all  transient  quantities,  the  free  surface  boundary  condi- 
tions are: 


P (r,()),0,u))  = P (r,cj)  , 

Z Z 0 

P,^(r,4),0,a3)  = Q (r,a))  , 

Z 0 


(3.1.5) 


where  P and  P are  the  Fourier  transformed  normal  and 
zz  zr 

shear  stress  components  in  a cylindrical  coordinate  system. 

The  P and  Q are  the  applied  stress  distributions. 

0 0 

The  elastic  waves  generated  by  loading  of  the  form 
(3.1.5)  on  the  free  surface  of  a multilayered  half space  can 
be  worked  out  analytically  using  techniques  similarly  applied 
by  Harkrider  [1964 J or  Bache  and  Harkrider  [1976]  for  the 
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burifed  source.  The  theory  is  presented  in  detail  in  our 
forthcoming  topical  report.  In  this  report  we  will  simply 
give  the  results. 

In  computing  surface  waves  for  the  ejecta  we  found 
the  maximum  amplitude  of  the  generated  surface  waves  to  be 
1-2  orders  of  magnitude  lower  than  for  the  surface  waves  from 
the  contained  portion  of  the  explosion.  The  amplitude  of 
the  displacement  spectrum  falls  off  rapidly  for  frequencies 
above  0.1  Hz  so  the  contribution  to  the  body  waves  is  quite 
a bit  less.  We  concluded,  therefore,  that  the  ejecta  con- 
tribution to  m^^  and  could  safely  be  ignored. 


3.2  BODY  WAVE  ATTENUATION  PROPERTIES  FOR  SELECTED  TRAVEL 


3.2.1  Introduction 

The  frequency  content  of  the  first  few  cycles  of  explo- 
sion P waves  recorded  at  teleseismic  distances  is  primarily 
controlled  by  the  shape  of  the  source  function,  the  delay 
time  between  P and  pP,  and  the  anelastic  properties  (t*  = T/Q, 
where  T is  the  travel  time  and  Q is  the  average  path  material 
quality  factor)  of  the  source-receiver  propagation  path.  The 
primary  objective  of  this  particular  study  is  to  determine 
variations  in  t*  between  U.S.  and  USSR  test  sites  and,  in 
particular,  the  effect  that  this  variation  produces  on  yield 
estimates  from  teleseismic  P waves.  As  a first  step  toward 
a solution  to  this  problem,  we  have  compared  the  frequency 
content  of  a large  population  of  NTS  and  Eurasian  explosions 
at  several  recording  stations  where  we  have  previous  experience 
in  predicting  explosion  P waves. 


3.2.2  Station  Locations  and  Experimental  Data 


Short-period  seismogreuns  of  explosion  P waves  were  ob- 
tained from  two  different  sources  for  this  experiment.  One 
source  was  the  Palmer  network  of  seismic  stations  located 
throughout  Alaska  and  the  Aleutian  Islands.  The  names,  codes 
and  locations  of  four  of  the  Palmer  stations  used  in  this 
experiment  are  given  in  Table  3.2.1. 


The  second  source  of  data  was  two  stations  of  the  Atomic 
Energy  Detection  System  (AEDS)  network,  hereafter  referred  to  as 
Sites  A and  B,  located  at  teleseismic  distances  from  the  NTS 
and  Eurasian  test  sites. 


Using  these  explosion  seismogreuns,  the  dominant  periods 
of  the  "b"  and  "d"  phases  of  the  short-period  P waves  were 
measured  for  each  event.  Examples  of  these  two  measurements, 
Tj^  and  T^,  are  shown  in  Figure  3.2.1. 
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As  can  be  seen  in  Figure  3.2.1,  is  based  on  a half- 
cycle  within  the  first  portion  of  the  P wave  while  is 
usually  based  on  a later  occurring  half-cycle  corresponding 
to  the  maximum  2unplitude  excursion  within  the  first  three  or 
four  cycles  of  the  P wave.  For  events  located  in  the  same 
source  region  and  recorded  at  a common  station  the  b phase 
has  been  found  to  be  a more  precise  measure  of  explosion 
coupling  than  the  d phase  [Bache,  et  al. , 1975]  . This  is 
attributed  to  the  fact  that  the  b portion  of  an  explosion  P 
wave  is  less  contaminated  by  later  arriving  phases  (i.e., 
pP,  mantle  arrivals). 


3.2.3  Results  of  Data  Analysis 

The  period  measurements,  Tj^  and  T^,  are  plotted  in 
Figures  3.2.2  through  3.2.15  versus  estimates  of  explosion 
yields  for  a large  population  of  NTS  and  Eurasian  events. 

The  yield  estimates  are  based  on  body  wave  magnitude  deter- 
minations and  are  taken  from  Klepinger  [1974]  and  Alewine, 
et  al . [1975]  . The  results  of  this  analysis  are  categorized 
and  presented  according  to  type  of  measurement  (Tj^  or  T^)  , 
recording  station  and  event  source  region. 

Some  comments  based  on  the  explosion  data  from  the 
Palmer  stations,  PMS  (Figures  3. 2. 2-3. 2. 5) , KDC  (Figures 
3. 2. 6-3. 2.7)  and  ADK  (Figures  3. 2. 8-3. 2. 9)  are  the  following: 

1.  The  dominant  periods  of  both  the  b and  d phases 
for  presumed  explosions  in  east  Kazakhstan  are 
shorter  than  the  corresponding  periods  for  NTS 
events  in  the  same  yield  range.  Mean  values  of 
Tj^  and  T^  are  shorter  by  approximately  0.5  sec 
at  PMS  and  KDC,  and  0.2  sec  at  ADK. 

2.  The  scatter  in  the  observed  Tj^  values  is  between 
a factor  of  1-1/2  and  2 less  than  the  scatter  in 
T^  values. 
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versus  Yield 


Tjj  versus  yield  for  large  population  of  NTS  and  Eurasian  (PNE's  and 
Novaya  Zemlya)  events  recorded  at  PMS. 


versus  Yield 


Figure  3.2.4.  T,  versus  yield  for  Seime  event  population  as  in  Figure 


versus  Yield 


same  event  population  as  in  Figure 


versus  Yield 


Station  of  the  Palmer  network 


versus  Yield 
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Figure  3.2.10.  Tt  versus  yield  for  NTS  and  Eurasian  (east  Kazakhstan)  events  recorded 
at  the  GIL  station  of  the  Palmer  network. 


versus  Yield 


versus  Yield 


Figure  3.2.12.  T,  versus  yield  for  NTS  (Yucca)  and  east  Kazakhstan  events  at  GIL 


PILEDRIVER 


Figure  3.2.13.  T versus  yield  for  explosions  recorded  at  SITE  A. 
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Figure  3.2.14.  T,  versus  yield  for  explosions  recorded  at  SITE 


Figure  3.2.15.  T.  versus  yield  for  explosions  recorded  at  SITE 
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3.  There  is  no  apparent  dependence  of  dominant  period, 
either  Tj^  or  T^,  on  event  yield  for  east 
Kazakhstan  explosions. 

4.  In  the  case  of  NTS  events,  there  is  no  observed 
dependence  of  period  on  yield  at  the  station  PMS. 

At  the  stations  KDC  eind  ADK,  however,  dominant 
periods  tend  to  increase  with  increasing  explosion 
yields. 

5.  At  PMS,  the  only  station  of  the  three  in  question 
for  which  data  from  Eurasian  explosions  not  lo- 
cated in  east  Kazakhstan  are  available,  there  is 
an  apparent  lengthening  of  dominant  period  with 
increasing  event  yield.  An  important  point  to 
note,  however,  is  that  these  events  occur  in  dif- 
ferent source  regions  at  significantly  different 
epicentral  distances  from  PMS  complicating  any 
interpretation. 

The  observations  from  the  remaining  Palmer  station, 

GIL,  are  included  here  to  point  up  the  problem  of  relying 
on  data  from  a single  station.  In  Figures  3.2.10-3.2.12, 
values  of  the  dominant  periods  are  plotted  versus  explosion 
yield  for  the  event  populations  considered  previously.  In 
contrast  to  the  pronounced  separation  of  event  populations 
seen  previously,  the  estimates  of  Tj^  and  T^  for  the  Eurasian 
(east  Kazakhstan  and  PNE's)  and  NTS  explosions  overlap 
significantly.  In  fact,  the  mean  Tj^  for  east  Kazakhstan 
events  is  slightly  longer  (0.85  sec)  than  that  for  NTS 
events  (0.75  sec).  For  PNE  versus  NTS  events,  the  mean 
values  of  T^^  are  approximately  equal.  The  apparent  reversal 
of  T^  for  east  Kazakhstan,  or  near  equality  for  the  PNE, 
versus  NTS  events  is  attributed  to  an  obvious  phasing  of  the 
GIL  seismograms  of  Eurasian  events.  This  phasing  (later 
arrival)  results  in  longer  apparent  periods  when  the  measure- 
ments are  made  according  to  the  convention  described  in 
Figure  3.2.1. 
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The  final  set  of  figures  shown  here  (Figures  3.2.13- 
3.2.15)  compare  dominant  periods  for  several  Eurasian  events 
with  a few  NTS  explosions  including,  in  particular,  the 
PILEDRIVER  event.  These  data  were  obtained  from  two  AEDS  stations 
at  teleseismic  distances  from  NTS.  Data  from  PILEDRIVER 
were  not  available  from  the  Palmer  stations  since  the  Palmer 
network  began  operating  during  September  1967,  more  than  one 
year  after  the  occurrence  of  PILEDRIVER. 

As  in  the  case  of  the  Palmer  stations  PMS,  KDC  and 
ADK,  the  estimates  of  Tj^  for  NTS  events,  including  PILEDRIVER, 
are  longer  than  Tj^  estimates  for  Eurasian  explosions  in  the 
same  yield  range.  A significant  point  about  the  PILEDRIVER 
explosion  is  the  fact  that  it  was  detonated  in  granite  and 
was  characterized  by  a fairly  short  pP-P  lag  time  compared 
to  other  NTS  events  of  comparable  yield.  It  is  clear  from 
the  results  in  Figures  3.2.13-3.2.15,  however,  that  this  ap- 
parently does  not  influence  the  observed  T.  estimates  rela- 

b 

tive  to  the  other  NTS  events. 

3.2.4  Conclusions 

As  noted  at  the  beginning  of  this  section  of  the  re- 
port, the  frequency  content  of  explosion  P waves  recorded  at 
teleseismic  distances  is  primarily  controlled  by  the  explo- 
sion source  function,  pP-P  lag  time  and  the  anelastic  proper- 
ties of  the  source-to-receiver  propagation  path.  Our  obser- 
vations indicate  that  P waves  from  Eurasian  events,  particularly 
those  originating  at  east  Kazakhstan,  are  of  significantly 
higher  frequency  than  P waves  from  NTS  events  recorded  at  the 
same  stations.  While  these  data  do  not  argue  conclusively 
for  differences  in  attenuation  beneath  the  NTS  and  Eurasian 
test  sites  being  the  cause  of  the  observed  period  (frequency) 
differences,  this  is  the  suggested  explanation  in  view  of 
the  large  sampling  of  different  event  depths  (pP-P  lag  times) 
and  source  media. 
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The  data  described  here  provide  important  constraints 
on  our  waveform  predictions  for  Eurasian  events.  Given  in- 
dependent estimates  of  pP-P  lag  times  (i.e.,  based  on  spectral 
nulls)  and  a range  of  predicted  source  functions,  we  will  be 
able  to  estimate  values  of  t*  for  the  different  test  site-to- 
receiver  paths  by  matching  the  observed  frequency  content 
of  seismograms. 
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DEVELOPMENT  OF  IMPROVED  TECHNIQUES  FOR  ONE-DIMENSIONAL 
NUMERICAL  SIMULATION  OF  EXPLOSIONS 


3.3.1  Introduction 

During  the  past  several  years  we  have  been  developing 
constitutive  models  and  stress  wave  propagation  techniques 
for  modeling  underground  nuclear  explosions  [e.g..  Cherry, 
et  al. , 1973;  1975].  The  computational  techniques  developed 
are  capable  of  simulating  the  explosion  induced  shock  wave 
from  the  near  source  high  pressure  regime  out  into  the  region 
where  the  material  response  is  linearly  elastic.  Most  of  our 
calculations  have  been  carried  out  in  a one-dimensional 
(spherically  symmetric)  geometry  for  which  we  use  the  computer 
program  we  call  SKIPPER. 

In  recent  months  we  have  needed  to  compute  the  equiva- 
lent elastic  source  for  explosions  in  highly  porous  materials 
exemplified  by  the  dry  tuffs  occurring  at  shallow  depths  at 
Yucca  Flat,  NTS.  Our  procedure  for  computing  the  equivalent 
elastic  source  is  to  monitor  the  displacement  field  in  the 
linear  elastic  field  and,  using  quasi-analytic  techniques, 
to  obtain  a reduced  displacement  potential  (RDP)  representa- 
tion of  the  field.  In  essence  we  invert  the  differential 
equation 


u(t-R/a)  = ^ 

r2  Ra 


(3.3.1) 


to  obtain  the  RDP,  4*. 

A satisfactory  RDP  must  exhibit  the  following  charac- 
teristics [Cherry,  et  al. , 1973]  : 

1.  The  RDP  should  be  invariant  with  distance,  R. 

2.  When  u « 0,  we  should  have  V « 0. 

3.  At  long  time  the  Y is  constant  ('F  'Vji  and 
u_  - R*^. 
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Typical  results  for  an  explosion  calculation  in  tuff 
with  14.5  percent  air-filled  voids  are  shown  in  Figure  3.3.1. 
Examining  the  figure,  it  is  clear  that  the  results  violate 
the  three  conditions  listed  above.  We  point  out  that  for 
materials  with  less  air-filled  porosity,  the  three  conditions 
are  satisfied  with  considerable  accuracy.  In  fact,  this  is 
required  to  match  analytic  solutions  such  as  that  by  Blake 
[1952] , a check  carried  out  at  the  time  SKIPPER  was  developed 
[Cherry,  et  al. , 1973] . 

In  view  of  the  failure  of  SKIPPER  to  properly  handle 
waves  in  highly  porous  materials,  the  constitutive  models 
were  reviewed  and  tests  against  analytical  results  were 
initiated  to  circumvent  this  difficulty.  The  results  of  this 
effort  are  described  in  the  sequel. 

3.3.2  An  Improved  Constitutive  Model  for  Wave  Propagation  In 

Porous  Materials 

A series  of  elastic  test  problems  were  initiated  to 
test  the  ability  of  SKIPPER  to  propagate  elastic  waves. 
Monitoring  the  pressure  versus  volume,  it  was  confirmed  that 
loading  and  unloading  occurred  on  the  same  path.  Further, 
the  numerical  results  agreed  very  well  with  analytical  solu- 
tions as  long  as  the  maximum  pressure  was  small  compared  to 
the  elastic  limit  pressure,  p^.  However,  the  differences 
became  noticeable  as  the  maximum  pressure  increased  to  a 
significant  fraction  of  p^. 

The  observations  stated  above  suggest  that  the  dif- 
ficulty is  in  the  fact  that  while  the  code  is  elastic,  it 
is  not  linear  elastic  as  the  p-V  relation  has  significant 
curvature  for  pressures  near  p^.  While  this  was  known,  the 
extent  of  the  effect  on  the  solution  was  not  known. 

The  problem  and  its  solution  are  summarized  graphically 
in  Figure  3.3.2.  In  (a)  we  show  the  p-V  curve  in  its  origi- 
nal form  for  porous  materials.  Even  though  the  material 
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loads  and  unloads  along  this  line  (elastic  behavior) , there 
is  significant  deviation  from  linearity  for  large  changes  in 
pressure.  For  our  elastic  test  problems  to  be  compared  to 
analytic  solutions,  we  merely  forced  the  p-V  relation  to  be 
linear  as  indicated  in  Figure  3.3.2b.  An  alternative  solu- 
tion is  to  choose  a lower  pressure  as  the  elastic  limit. 
Although  this  method  would  clearly  work,  it  would  increase 
the  cost  of  the  calculations  and  probably  would  not  pro- 
duce results  significantly  different  from  the  solution 
chosen. 

In  a full  explosion  calculation  it  is  necessary  to 
have  a p-V  relationship  that  extends  into  the  nonlinear  re- 
gime. The  adjusted  p-V  relation  for  the  porous  tuff  calcula- 
tion of  Figure  3.3.1  is  shown  in  Figure  3.3.2c.  The  data 

were  made  linear  to  p and  a smooth  transition  between  p 

*^e  e 

and  the  nonlinear  portion  of  the  crust  curve  was  introduced. 
This  modification  in  the  material  behavior  does  not  contra- 
dict available  laboratory  data. 

The  calculation  described  in  Figure  3.3.1  was  repeated 
with  the  minor  modification  described  above.  The  results 
are  shown  in  Figure  3.3.3.  We  see  that  now  the  calculations 
satisfy  the  three  conditions  of  Section  3.3.1  quite  well. 

Our  ultimate  objective  in  computing  RDF's  for  under- 
ground explosions  is  to  determine  the  effect  of  local  material 
properties  on  such  gross  measures  of  seismic  energy  as  m^ 
and  Mg.  The  question  is,  how  much  has  our  adjustment  of  the 
constitutive  model  affected  these  quantities.  The  answer  is 
apparent  from  Figure  3.3.4  where  we  show  the  aunplitude  of  the 

1 ^ I 

Fourier  transformed  reduced  velocity  potential  (|4'(a))|)  for 
the  two  calculations.  The  source  functions  are  appropriate 
for  a yield  of  20  kt. 

The  nij^  is  controlled  by  the  portion  of  the  spectrum 

from  t 0.8  - 2.0  Hz.  The  M is  sensitive  to  the  long  period 

s 

spectrum  with  frequencies  from  t 0.04  - 0.1  seconds.  Com- 
paring the  two  source  functions  we  see  that  in  this  case 
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relatively  minor  chcmges  result  from  the  introduction  of  the 
altered  p-V  relation.  We  expect  the  nij^  and  to  increase 
slightly,  no  more  than  0.1  magnitude  units,  with  the  improved 
constitutive  model.  However,  the  RDP  in  the  first  calculation 
is  range  dependent  and  much  greater  differences  could  result 
if  we  compared  to  the  RDP  from  another  range. 


In  summary,  we  have  found  that  deviations  from  linear 
elasticity  can  significantly  compromise  the  ability  to 
compute  a satisfactory  RDP  for  1-D  explosion  calculations. 

The  problem  can  be  hidden  in  the  sense  that  it  only  becomes 
important  for  certain  material  properties,  in  this  case  high 
porosity.  Only  minor  modifications  to  the  constitutive 
model  were  required  to  fix  the  problem.  The  use  of  laboratory 
data  to  select  the  rock  properties  input  to  the  program  re- 
mains unchanged.  Finally,  the  SKIPPER  code  was  once  again 
tested  against  analytical  solutions  in  the  elastic  regime,  a 
test  that  should  be  periodically  done  for  programs  being 
modified  and  improved. 
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3.4  THE  DEPENDENCE  OF  BODY  WAVE  MAGNITUDE  ON  YIELD  FOR 

UNDERGROUND  EXPLOSIONS  IN  SALT 

A theoretical  study  was  conducted  in  which  the  objec- 
tive was  to  compare  explosions  in  salt  in  Eurasia  to  similar, 
though  hypothetical,  events  at  NTS  and  to  granite  events  at 
NTS.  The  results  of  this  study  were  reported  in  detail  in  a 
topical  report  by  Bache,  Cherry  and  Mason  [1976].  In  this 
report  we  reproduce  the  main  results  from  the  topical  report 
and  give  some  additional  results:  theoretical  spectra  for 
several  of  the  calculations. 

The  equivalent  elastic  source,  the  reduced  displacement 
potential,  was  computed  for  four  explosions  in  salt  in  which 
only  the  overburden  pressure  was  varied.  The  source  spectra 
are  shown  in  Figure  3.4.1.  The  relationship  between  the  'i' (to) 
and  the  far-field  displacement  spectrum  is  given  by  Eq. 

(3.1.2) . 

From  the  spectra  of  Figure  3.4.1,  we  see  that  the  long 
period  level  is  mildly  sensitive  to  depth,  varying  by  only 
20  percent  over  the  entire  range.  However,  the  short  period 
level  varies  by  a factor  of  2.2  from  largest  to  smallest. 

This  is  because  of  the  marked  peaking  of  the  source  function 
at  shallow  depths.  This  peaking  is  due  to  the  opening  of 
tensile  cracks  in  the  material  at  low  confining  pressures. 

The  sources  of  Figure  3.4.1  were  used  to  compute  syn- 
thetic short  period  seismograms  at  teleseismic  distances. 

A key  parameter  for  these  calculations  is  the  attenuation 
parameter  t*.  For  paths  such  as  NTS-Alaska  that  we  have 
studied  extensively,  we  believe  t*  = 1.0.  Because  we  believe 
t*  is  lower  for  Eurasian  travel  paths  we  rather  arbitrarily 
have  chosen  t*  = 0.6  for  these  calculations. 

The  computed  synthetic  seismograms  for  a particular 
range  are  shown  in  Figure  3.4.2.  The  source  function,  yield 
and  burial  depth  for  each  calculation  are  printed  on  the 
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Frequency  (Hz) 

Figure  3.4.1.  Source  functions  for  salt  with  variable  over- 
burden (indicated  by  the  depth  H,  in  km).  Thi 
frequency  axis  is  scaled  to  160  kt»  the  ampli' 
tude  axis  to  0.02  kt. 
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records.  Also  appearing  on  the  records  are  the  values  and 
the  period  of  the  cycle  at  which  is  measured. 

In  Figure  3.4.3  we  plot  nij^  versus  log  yield  (W)  from 
the  synthetic  seismogreuns.  The  best  least  squares  linear  fit 
to  the  data  was  computed  for  the  total  group  and  two  sub- 
groups depending  on  the  scaled  depth  of  burial. 

We  next  change  the  path  parameters  to  be  appropriate  to 
the  NTS-Alaska  travel  path  that  we  have  previously  calibrated. 
The  most  important  change  is  now  t*  = 1.06.  Synthetic  seismo- 
greims  were  computed  for  five  of  the  twelve  events  of  Figure 
3.4.2.  The  seismogrcims  are  shown  in  Figure  3.4.4. 

We  are  also  interested  in  comparing  lOj^-yield  for  the 
salt  events  to  that  for  explosions  in  granite.  We  have  pre- 
viously [Bache,  et  al . , 1975]  computed  a theoretical  m^^  for 
PILEDRIVER  that  agrees  quite  well  with  observations.  We 
expect  the  nvj^-yield  curve  for  NTS  granite  explosions  would 
be  close  to  a curve  of  unit  slope  through  the  PILEDRIVER 
datum. 

The  m^^-yield  values  for  salt-Eurasian  travel  path, 
salt-NTS  travel  path  and  PILDRIVER  granite-NTS  travel  path 
are  compared  in  Figure  3.4.5.  Three  lines  of  unit  slope  are 
shown  through  the  data.  The  uppermost  line  is  the  least 
squares  fit  to  the  salt-Eurasian  path  events  of  Figure  3.1.2. 
The  lowest  line  is  0.55  magnitude  units  below  the  first. 

This  is  representative  of  the  difference  between  the  salt- 
Eurasian  path  and  salt-NTS  path  mj^  values  if  the  t*  change 
from  0.6  to  1.06  is  approximately  correct. 

The  unit  slope  line  through  the  PILEDRIVER  theoretical 
m^  falls  almost  equidistant  between  the  lines  for  the  two 
populations  of  salt  events.  That  is,  events  like  PILEDRIVER 
are  expected  to  give  m|^  values  that  are,  on  the  average, 

0.1  - 0.3  magnitude  units  above  those  for  similar  NTS  events 
in  salt.  On  the  other  hand,  the  granite  events  give  nVj^  values 
that  are  about  0.1  - 0.3  units  lower  than  Eurasian  salt  events. 


Figure  3.4.4,  Theoretical  seismogreuns  for  hypothetical  events 
in  salt  at  NTS.  The  source  parameters  are  the 
same  as  for  the  corresponding  Eurasian  salt 
events  studied  in  previous  sections.  The  path 
par^uneters  are  appropriate  for  an  NTS-Alaska 
path. 
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Figure  3.4.5.  Comparison  of  mt-yield  for  salt  events  in  Eurasia  and  at  NTS  to  granite 
explosions  at  NTS. 
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Not  presented  in  the  Bache,  Cherry  and  Mason  [1976] 
topical  report  but  of  some  interest  are  typical  displacement 
spectra  for  our  theoretical  salt  events.  We  select  four  of 
the  events  for  which  seismograms  were  presented  in  Figure 
3.4.2.  Displacement  spectra  for  these  calculations  are  shown 
in  Figures  3.4.6  - 3.4.7.  In  Figure  3.4.6a  we  show  the  P 
wave  velocity-depth  model  used  for  the  crust  in  the  source 
region.  The  spectra  in  the  figure  are  for  the  wave  as  it 
exits  the  bottom  of  the  crust  at  the  source.  The  highest 
frequency  for  the  spectral  computations  was  5 Hz. 

The  most  important  feature  controlling  the  shape  of 
the  spectra  in  Figure  3.4.6  is  the  free  surface  reflection 
or  pP  phase.  This  point  is  illustrated  in  Figure  3.4.5b 
by  the  spectrum  for  one  of  the  events  with  the  crustal  model 
modified  to  allow  the  salt  layer  to  extend  upwards  to  infinity. 
Then  the  only  energy  in  the  spectrum  is  that  which  leaves 
the  source  downward . 


In  Figure  3.4.7  we  show  the  spectra  at  the  receiver, 
with  and  without  the  instrument,  for  the  same  four  events. 
Once  again  given  for  comparison  are  the  spectra  for  the  40  kt 
event  with  the  layering  above  the  source  removed.  For  these 
spectral  calculations  the  Nyquist  frequency  is  at  about  0.1 
Hz  and  the  portion  of  the  plot  below  this  frequency  should  be 


Salt 


Figure  3.4.6.  Displacement  spectra  at  the  base  of  the  crust  below  the  source  fo 
four  of  the  salt  events  for  which  seismograms  are  given  in  Figure 
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Figure  3.4.7,  Ground  motion  displacement  spectra  for  the  events 
of  Figure  3,4.6  with  and  without  the  instrument 
response  included. 
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3.5  A LABORATORY  SIMULATION  OF  UNDERGROUND  EXPLOSIONS 

Carried  over  from  our  previous  contract  (Project  VT/ 
6712)  is  a progreun  in  which  small-scale  modeling  experiments 
of  underground  explosions  are  being  carried  out  and  the  re- 
sults compared  to  computer  simulations  of  the  experiments. 
During  the  past  quarter,  several  successful  model  experiments 
have  been  conducted  and  some  results  are  presented  here. 

The  objective  of  the  experiments  is  to  obtain  displace- 
ment time  histories  on  the  free  surface  of  a uniform  medium 
of  well-known  properties.  The  measurements  are  made  on  the 
free  surface  at  large  enough  distances  to  ensure  that  the 
response  is  linear  elastic. 

The  test  medium  was  a "high-early"  concrete  of  density 
2.52  g/cm*,  unconfined  compressive  strength  3.45  x 10 ’ Pa 
(345  bar) , low  porosity  and  high  water  saturation.  The  ob- 
served sound  speed  was  3.4  x 10®  cm/s.  Each  explosive  charge* 
consisted  of  1/4  g of  PETN,  packed  to  a density  of  1.0  g/cm* 
inside  a lucite  shell  of  radius  0.56  cm.  Great  care  was 
exercised  in  the  construction  of  each  charge  to  ensure  that 
the  resulting  detonation  was  spherically  symmetric. 

The  principal  sensors  were  eddy-current  displacement 
gauges.**  Each  unit's  output  varies  linearly  with  distance 
between  the  probe's  tip  and  a small  piece  of  aluminum  foil 
fixed  with  epoxy  to  the  concrete  surface.  These  transducers 
have  rather  low  output  levels,  2 mV/ym,  but  they  have  the 
advantages  that  they  do  not  touch  the  concrete  and  their 
response  extends  to  50  kHz  (3  db) . To  supplement  the  dis- 
placement gauges,  we  also  used  two  accelerometers  with  ranges 
of  + 50,000  g's  (5  x lo*  m/s*),  sensitivities  of  10  VV/g  and 


* Supplied  by  Stanford  Research  Institute,  Menlo  Park,  CA. 

Supplied  by  Kaman  Sciences,  Colorado  Springs,  Colorado. 

^ Supplied  by  Endevco,  San  Juan  Capistrano,  CA. 
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resonant  frequencies  above  150  kHz.  The  low  mass  of  the 
accelerometers  allowed  them  to  be  epoxied  to  the  concrete 
directly. 

We  prepare  the  concrete  in  the  form  of  right  circular 
cylinders  of  diameter  120  cm  and  thickness  60  cm.  The  motion 
transducers  were  arrayed  over  one  end  of  the  cylinders  and 
the  charges  were  located  31  cm  below  those  surfaces. 

In  one  experiment,  two  displacement  gauges  and  an 
accelerometer  were  in  close  proximity  and  directly  over  the 
charge.  Our  objective  was  to  see  if  the  three  independent 
transducers  would  give  similar  responses  to  the  seismic  motion. 
In  fact,  the  responses  were  quite  similar,  showing  an  arrival 
time  of  92  ys,  risetime  (in  displacement)  of  15  ys,  peak  mo- 
tion of  5 ym  and  full  width  at  half  maximum  of  about  20  ys. 

For  the  subsequent  tests,  the  gauges  were  located  at 
the  vertices  of  a one  foot  equilateral  triangle;  horizontal 
(toward  vertex  #2)  and  vertical  displacement  at  vertex  #1; 
vertical  displacement  and  acceleration  at  vertex  #2;  and  ver- 
tical acceleration  at  vertex  #3.  Charges  were  placed  31  cm 
below  each  vertex  and  fired  one  at  a time.  Some  typical  re- 
sults are  shown  in  Figures  3.5.1  and  3.5.2.*  While  there  are 
several  reasons  to  believe  that  the  double  integral  of  an 
accleration  signal  is  not  equal  to  the  true  displacement, 
especially  at  late  times,  the  acceleration  data  do  generally 
confirm  the  shape  and  amplitude  of  the  main  displacement  peak. 

As  might  be  expected.  Figure  3.5.1  shows  that  no  hori- 
'zontal  motion  (HI)  occurs  until  the  arrival  of  the  more  slowly 
moving  shear  wave  at  about  160  ys.  All  of  the  data  records 
also  Indicate  that  after  the  main  displacement  peak,  the  sur- 
face continues  to  have  a small  but  consistently  positive 


The  error  bars  shown  include  the  uncertainties  due  to  instru- 
ment and  recorder  calibrations  and  electrical  noise. 
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Data  like  that  in  the  previous  figure  except 
that  now  the  charge  is  located  under  vertex  2 
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displacement  >5ver  an  interval  of  at  least  100  ys;  see  VI  in 
Figure  3.5.1  and  V2  in  Figure  3.5.2. 

The  experiments  described  above  will  be  modeled  with 
our  finite  difference  programs  and  will  thus  provide  an  ex- 
cellent normalization  of  the  theoretical  techniques.  Further 
experiments  are  also  being  performed  in  which  the  charge  em- 
placement will  be  shallow  and  cratering  will  result. 
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